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We investigate the pressure-induced metal-insulator transition from diamond to /3-tin in bulk Sili- 
con, using quantum Monte Carlo (QMC) and density functional theory (DFT) approaches. We show 
that it is possible to efficiently describe many-body effects, using a variational wave function with 
an optimized Jastrow factor and a Slater determinant. Variational results are obtained with a small 
computational cost and are further improved by performing diffusion Monte Carlo calculations and 
an explicit optimization of molecular orbitals in the determinant. Finite temperature corrections 
and zero point motion effects are included by calculating phonon dispersions in both phases at the 
DFT level. Our results indicate that the theoretical QMC (DFT) transition pressure is significantly 
larger (smaller) than the accepted experimental value. We discuss the limitation of DFT approaches 
due to the choice of the exchange and correlation functionals and the difficulty to determine con- 
sistent pseudopotentials within the QMC framework, a limitation that may significantly affect the 
accuracy of the technique. 

PACS numbers: 



I. INTRODUCTION 

The prediction of material behavior under pressure is 
of great relevance in several branches of science, ranging 
form material science to planetary physicsA. The exper- 
imental determination of pressure effects on the phase 
diagram can be rather complicated even for simple inor- 
ganic crystal, due to the related extreme conditions or the 
existence of subtle phenomena, such as hysteresis. In this 
respect, ab-initio calculations are complementary meth- 
ods for determining phase diagrams and understanding 
material phases in a large range of pressures and temper- 
atures. 

The density functional theory (DFT) has been widely 
used for describing material behavior under pressure. 
However, when a transition is accompanied by a dras- 
tic change in the electronic structure, non-canceling er- 
rors in the two phases can lead to a significant bias 
in the predicted transition pressure. That is the case 
of bulk Silicon (Si), where the diamond-to- /3-tin tran- 
sition is associated with a semiconductor-to-semimetal 
electronic change. For this reason, the Si diamond-to-/3- 
tin transition has been used for testing and benchmarking 
new ab-initio numerical approaches for extended systems. 
Its first order nature makes it a difficult case for the 
experiments 2 " 4 . The accepted experimental values for 
the transition pressure are in between 10 and 12.5 GPa 
at room temperature, the difference could be ascribed to 
non hydrostatic conditions^, non quenched metastable 
phases^, and presence of lattice defects 7 -. 

Theoretical calculations based on DFT strongly de- 
pend on the choice of exchange and correlation (XC) 
functional. Calculations based on the local density ap- 



proximation (LDA) predict a zero temperature transition 
pressure pt in the range of 7.3^-8.3^ GPa (this difference 
could be explained by the type of pseudopotential (PP) 
used). The generalized gradient approximation (GGA) 
leads to pt in the range of 12.3i^-13.5^ GPa with the 
Perdew-Wang functional 11 . Moreover, in the case of the 
PBE functional 12 , a value of 10.2 GPa is obtained in 
Ref. [^. Since the PBE functional fulfills a number of ex- 
act DFT properties, this should be considered the state- 
of-the-art among the most accurate ab-initio functionals. 
A careful investigation of the effect of the XC functional 
on the transition pressure was done recently in Ref. Il3l . 
where the authors show that the inclusion of non-local 
exchange in the XC functional leads to a significant im- 
provement of the estimate of the transition pressure. All 
the calculations are performed at zero temperature and 
therefore a comparison with experiments is possible only 
after including zero-point motion and finite temperature 
effects. This accounts for a significant reduction of the 
transition pressure, as estimated in Refill by means of 
the LDA functional, with a correction larger than lGpa. 
Moreover the explicit inclusion of non linear core cor- 
rections (NLCC) in pseudopotential calculations further 
reduces the transition pressure. 

Quantum Monte Carlo (QMC) methods can be an 
alternative to DFT-based approaches. In the past 
years, many authors have shown practical applications 
of QMC methods for computing the energetics of ex- 
tended systems 15, — , and predicting crystal phases under 
pressur o 17 ! 18 . In a early work, Alfe et al~ used diffusion 
Monte Carlo (DMC) for investigating the Si diamond- 
to- /3-tin transition. They calculated a QMC transition 
pressure of 16.5 GPa, namely 4-5 GPa larger than the ex- 
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perimental range. The discrepancy was attributed to the 
fixed-node (FN) approximation 1 ^, since the other source 
of errors (time step error, pseudopotential locality error, 
size effects) were considered negligible. More recently, 
Purwanto et al. 2 ® obtained a transition pressure of 12.6 
using the auxiliary-field QMC (AFQMC) method with 
the phaseless approximation to cure the sign problem, 
whose bias has been shown to be very small^i. The very 
recent DMC calculation in Ref. [HI, carried out with the 
most advanced optimization 2 ^ and size-extrapolations 3 - 
methods, gives a transition pressure of 14 GPa in sub- 
stantial agreement with the AFQMC one, although the 
latest DMC value is affected by a larger uncertainty 
(lGPa). 

In the present work, we address the problem of the 
Si diamond-to- /3-tin zero temperature transition by per- 
forming variational Monte Carlo (VMC) calculations for 
the total energy of the two crystal structures. VMC is 
not affected by the sign problem and has proven to be a 
reliable approach for several system s 15 ' 16 ' 24 ' 25 . The ac- 
curacy of the method depends entirely on the choice of 
the variational wave function, and the capability of find- 
ing its optimal form. In this work, we show that with 
a relative simple parameterization of the wave function 
(WF) we are able to describe correlation effects across 
the metal-insulator transition. Our wave function is a 
product of a Jastrow factor and a Slater determinant. 
The variational optimization of the Jastrow factor has a 
relative small computational cost^ 2 -, and it reveals an ac- 
curate way to build-up correlation effects starting from 
an LDA calculation. 

The paper is organized as follows. In Sec. HI] we re- 
view the properties of our variational wave function, we 
present a systematic study of the basis set used for the 
calculations, and we discuss the correction of the finite 
size errors. In Sec. IIIII we review sources of errors not 
directly estimated at the QMC level, such as the quan- 
tum and thermal lattice energies, and the accuracy of the 
PP's. This will help us to make a fair comparison be- 
tween our results and the most significant experimental 
and theoretical findings reported in literature. In Sec. lIVI 
we report our results, while in Sec. [V]we draw our con- 
clusions. 



II. COMPUTATIONAL DETAILS 

In our investigation we have carried out DFT calcu- 
lations with the Quantum Espresso (QE)2£, Wien2K 27 
and QbooP 8 packages, and we used the TurboRVB code 29 
to perform variational Monte Carlo (VMC) and Lattice 
Regularized Diffusion Monte Carlo (LRDMC) 3 ^ calcula- 
tions. In both DFT and QMC the system is described 
by an effective Hamiltonian in the Born-Oppcnhcimer 
approximation (without quantum effects of lattice vibra- 
tions) with the core electrons treated neon-core PP. 

For extensive volume (pressure) dependent calcula- 
tions, we used the relativistic Hartree-Fock PP's gen- 



erated by Trail and Needs^I, employed also in pre- 
vious QMC calculations of Refs. (fl3||l8[). To check 
the impact of the PP approximation we carried out 
single-volume calculations using the Burkatzki-Filippi- 
Dolg energy-adjusted PP's 3 ^ and the ones generated by 
the atomic PWSCF code 2 ^ with the Troullier-Martins 
construction 3 ^ from relativistic LDA and PBE atomic 
calculations. The PP's are expanded in a semilocal form 
in terms of Gaussians^ 6 -. 



A. Wave function and localized basis set 



Our many-body wave function is the product of a 
Slater determinant and a many-body Jastrow factor. Al- 
ternatively, in the case of open shell systems, we use 
the Jastrow-correlated Antisymmetrized Geminal Power 
(JAGP)2£. 

The electron correlation is included in our wave 
function through the Jastrow factor J(ri,-- - ,rjv) = 
]~[ exp(/(ri, Tj)), where /(r, r') is assumed to depend 

i<j 

only upon two-electron coordinates, and N is the total 
number of electrons. The function / is expanded in a 
basis of Gaussian atomic orbitals fa, such that it reads: 



/(r,rO = 5>,<Mr)^(r'). (1) 



The convergence of this expansion is improved by adding 
an homogeneous term and a one-body contribution, thus 
satisfying the electron-electron and the electron-ion cusp 
conditions, respectivel y 25 ' 34 . The basis set used for the 
Jastrow includes 2s2p Gaussian orbitals. 

For system with large number of electrons, we improve 
the efficiency of the optimization procedure and its com- 
putational cost, adopting an explicit parameterization of 
the Jastrow factor. Given two generic orbitals fa and fa, 
the variational coefficient g^ = ffy(Rij) in Eq.Q] depends 
only on the orbital symmetry (in this case either s or 
p) and the distance vector Ry between the correspond- 
ing atomic centers, gij is optimized without constraint 
when the two orbitals are localized on the same atom (i.e. 
Ry = 0). Otherwise gij (Rj j) is parameterized in a way 
to recover an isotropic large distance correlation. For 
the sake of clarity, we define as the vector containing 
the three x, y, z p-orbital components centered at a given 
atomic position R J; and with Sj we indicate the s-wave 
orbital located at the same position. By this notation we 
can write four possible isotropic invariant contributions 



3 



for Ry 7^ 0, such that: 

l— Si ,p ; 
m=Sj ,pj 

+ f3(Rij){rij ■ 4> Pi )(rij ■ 4> Pj ) 

(2) 

where Rij = |Rjj|, and r,y = Ry /i?y is the unit vector 
connecting the atomic centers i and j. The functions 
fp(Rij) are polynomials which read: 

/ p (i? 2 ,) = CI log(i^) + J] C£R£ fc for p = 1, 2 

fc=i 

fc=3 

/p(%) - forp = 3,4 (3) 

fe=i 

The scalar functions in Eq. [3] depend only on 12 varia- 
tional parameters, optimized via energy minimization 35 . 

We verified the validity of the chosen parameterization 
(at long and short distances), by a direct comparison 
with the case of a fully optimized Jastrow factor (i.e. 
without parameterization) . The expression in Eq. [5] can 
be appropriate for physical long distance behaviors of 
the Jastrow factor, including the one recently speculated 
for describing the Mott insulator, and containing a long 
range \og(Rij) term 3 ^. 

The Slater determinant is obtained with N/2 molecu- 
lar orbitals ipj (r) , each doubly occupied by opposite spin 
electrons. The orbitals ipj( r ) are expanded in a Gaussian 
single-particle basis set {4>i\, centered on atomic nuclei, 
i.e. ipj{r) = J2i Aji0i(r). The Slater determinant is build 
from LDA orbitals. LDA calculations are performed with 
a periodic Gaussian basis set (see Appendix [S] for its 
definition) by using the DFT code included in the Tur- 
boRVB package 29 . This allows us to perform an efficient 
DFT calculation in exactly the same basis used in QMC 
and without employing the so called Kleinman-Bylander 
approximation on the PP's. 

We carefully tested the effects of the basis set extrap- 
olation on the total energy for both diamond and /3- 
tin geometries. Following the systematically convergent 
method for accurate total energy calculations recently in- 
troduced in Ref.(H3), we have used a tempered basis set, 
where the Gaussian exponents Z± are defined as Z$ = a/3 t 
for % = 0, . . . , n — 1 with a = Z mul and a/3™ = Z max . 
The parameters n, Z min and Z max are free. We verified 
that our basis set parameterization guarantees the same 
accuracy in both metallic and insulating phases for all 
investigated pressures. 

For a given angular momentum I (< 4), we fixed the 
maximum number of Gaussians according to the formula 
rii = no — 21, inspired by the correlated consistent basis 
set approac h 38 ! 39 . It follows that the maximum number 



n of exponents is used only on the s-wave channel, where 
rio = n. We studied the basis set extrapolation with 
respect to n by fixing Z min = 0.05 and Z max = 10. 
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FIG. 1: Accuracy of the DFT energy as a function of the ba- 
sis set extension n(= no) defined in the text. The number of 
Gaussians in each angular channel I is given by m — no — 21, 
as explained in the text. In the plot, we indicate with an 
arrow the value of n when the corresponding angular chan- 
nel has been switched on in the basis set expansion. The 
red open squares refer to the metallic /3-tin phase at volume 
13. 081 A 3 /Si and c/a — 0.54, while the blue triangles to the 
diamond insulating phase at volume 20.036A 3 /5'i. 

In Fig. [TJ we show the convergence of the DFT energy 
as a function of n{= hq) for a system of 8 Si atoms. 
We report the accuracy at 13.08lA 3 /5i and c/a — 0.54 
for the /3-tin phase, and at volume 20.036A 3 /S'i for the 
diamond phase. The accuracy in the energy is estimated 
by using a fully converged reference energy of a plane- 
wave calculation with 100 Ry kinetic energy cutoff. This 
was obtained with the Qbox package 2 ^, by using the same 
PP's in the semilocal form as the ones used in our QMC 
calculations. An accuracy of 0.01 eV/Si is sufficient to 
determine the equation of state (EOS) with an error well 
below the experimental uncertainty. From this analysis 
it turns out that n > 7 and the inclusion of d-orbitals 
in the basis set guarantee an accuracy of 0.007 eV/Si on 
the energy difference between the two phases. 

In the following, VMC and LRDMC production runs 
are performed using a basis set with n = 12, ni = 6, and 
n% — 4. The basis set exponents have been optimized at 
the DFT level. We found the optimal parameters Z m i n — 
0.05 and Z max = 3.25, that minimize the DFT energy in 
both the diamond and j3— tin phases. 

All Jastrow parameters, including exponents, are ob- 
tained by means of VMC energy minimization 35 . Pos- 
sible correlation effects not included by our Jastrow pa- 
rameterization are recovered by performing LRDMC cal- 
culations. As a projector method, LRDMC allows one to 
obtain the best variational wave function with the same 
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nodal structure as the initial variational wave function, 
the so called fixed node approximation (FN), giving an 
upper bound of the true ground state energy even with 
non local pseudopotential a 30 ' 40 . 



B. Finite-size errors 

Contrary to standard DFT methods, QMC calcula- 
tions have to be performed on a supercell. Therefore, 
finite size (FS) effects can be a relevant source of error in 
QMC calculations. Several methods have been proposed 
for correcting FS errors. One source of FS errors arises 
from the kinetic and Hartree term and can be treated by 
standard DFT approach with fc— point sampling. This is 
a genuine one-body contribution, and can be corrected by 
#dft _ E% FT , i.e. the difference between the DFT en- 
ergy per atom with a fully converged k— point mesh and 
the energy per atom of the supercell with finite volume 
and number of electrons N. 

The other source of errors (two-body terms) is related 
to the finite size effects of the exchange and correlation 
(XC) functionals, not explicitly included in E® FT . We 
calculate the two-body term corrections using the func- 
tional proposed by Kwee, Zhang and Krakauer (KZK) 23 . 
In Ref. |23| the authors proposed to estimate this type 
of FS error within the LDA framework, where the ex- 
change and correlation energy functional is replaced by 
the LDA functional parameterized for a finite system, 
which keeps an explicit dependence on the number of 
particles. Therefore, the total one- and two-body correc- 
tion is given by A KZK = £ DFT - i^ F K T ZK , where E^ ZK 
is the DFT energy computed with the KZK functional 
for N electrons. 

We observe that FS errors can be particularly relevant 
for open shell metallic systems. The AGP wave function 
approach allows to include many determinants in a ef- 
fective way, removing the degeneracy of the open shell 
by an appropriate fractional occupation of the degener- 
ate levels. Within LDA, by using a negligible smearing 
in the occupation of the KS energy levels, the degener- 
ate orbitals containing the HOMO are partially occupied 
with the same charge, and can be used consistently in 
the AGP wave function. 

For alleviating the effects of the one-body terms, we 
performed VMC and LRDMC calculations averaging the 
KZK corrected energies over the two most symmetric 
points (r and M). Previous DMC calculation s 13 ' 18 were 
performed only at the T point, leading to more pro- 
nounced size effects. In Tab. U we report the energy for 
different sizes together with KZK corrections, while in 
Fig. [5] we show how important is to average over the 
T (PBC) and M (APBC) points to reduce considerably 
the error in the finite size extrapolation 5 -!. By averaging 
over PBC and APBC boundary conditions we reach an 
accuracy of 5 meV/atom with 64 atoms, well below the 
magnitude of other systematic errors, as we will see in 
Sec. IIIII Therefore, at variance with Ref. [HI, where a 
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FIG. 2: Finite size extrapolation for the /3-tin (V = 15 A 3 /Si, 
c/a=0.55) and the diamond (V = 19.949 A 3 /Si) phases in 
the upper and lower panel, respectively. The data points cor- 
respond to the values reported in Tab. [I] obtained with the 
Trail-Needs pseudopotentials. In the plot the energy per atom 
is shown as a function of 1 /N at oms , where iV a toms is the num- 
ber of Si in the crystal supercell. Note the slopes for a given 
boundary condition (either PBC or APBC) in the two cases 
differ by a factor of 10, making the extrapolation much harder 
for the metallic phase. 

single k-point was adopted in the size extrapolations, the 
finite size bias is not the largest error in our calculations. 

III. IMPACT OF VARIOUS APPROXIMATIONS 
A. Pseudopotential approximation 

The change in Silicon coordination number (from 4 to 
6), related to the structural transition from diamond to 
/3-tin geometries, may affect the transferability of Si PP's 
in the two phases. 

We verified the impact of PP's on our final results, by 
estimating the transition pressure at the DFT-LDA level 
with different norm conserving PP's and with the projec- 
tor augmented-wave (PAW) metho d 41 ' 42 ! 58 The results 
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N atoms 


pgMU a KZK 
-^PBC + A PBC 


a KZK 
^PBC 


ttiC^IVKJ a KZK 
-^APBC + A APBC 


a KZK 
^APBC 


. j-jUiviC , t-iOMU , a KZK , a KZK \ 
l^PBC + -^APBC l A PBC + ^APBCj/^ 


/ a KZK a KZK \ /o 
l^PBC + ^APBCj/2 


64 ,8— tin 


-106.3063(17) 


-.0654 


-106.3705(17) 


.2214 


-106.3384(12) 


.07799 


96 /3-tin 


-106.3103(17) 


-0.0584 


-106.3572(20) 


0.0865 


-106.3338(13) 


0.01409 


256 /3-tin a 


-106.3395(23) 


-.0273 


-106.3417(32) 


.08855 


-106.3406(22) 


.03064 


64 diamond 


-106.7871(23) 


-0.0329 


-106.7791(23) 


-0.0044 


-106.7831(16) 


-0.0187 


8 diamond 


-106.6864(14) 


-.6700 


-106.6786(14) 


-.6632 


-106.6825(10) 


-0.6666 



TABLE I: Energy per atom (in eV) for various system sizes in the metallic /3— tin (V = 15 A 3 /Si, c/a=0.55) and insulating 
diamond (V = 19.949 A 3 /Si) phases with the Trail-Needs pseudopotentials. £-S?p BC is the variational QMC energy with the 

given PBC or APBC boundary conditions. Aou PBG is the total one-body and two-body correction computed by means of the 
KZK energy functional, with the same boundary conditions. 

a Corrected by — 0.0082(13)eV/atom to take into account that in this case for the long distance tail of the Jastrow a few 
variational parameters ansatz was adopted, as described in Sec. Ill Al The same form was used in the 64— Si case to 
estimate this correction. 



have been compared with all electron calculations ob- 
tained with the Wien2k2£ code. In Fig. [31 we report the 
EOS obtained with norm conserving PP's generated by 
the Troullier-Martins method with scalar relativistic cor- 
rections and different cutoff radii r c , the data obtained 
by the PAW method, and the reference all-electron re- 
sults. The PP DFT calculations have been done with 
the PWSCF code2&. We worked with a plane wave cut- 
off of 50 Ry and a charge density cutoff of 200 Ry. The 
number of non-equivalent k-points in the Brillouin zone 
is 160 for the /3-tin phase, and 80 for the diamond struc- 
ture. We checked that those parameters give converged 
DFT results with a Gaussian broadening of 0.01 Ry of the 
Fermi surfaced On the other hand, the Wien2k calcula- 
tions are PP error free, and therefore they can be used to 
check the PP accuracy. They have been performed with 
an equivalent Brillouin zone integration over a 17x 17x 17 
k-point mesh, a muffin-tin radius Rmt — 1-90 an (an is 
the Bohr radius), and a plane wave cutoff fc max given by 
^MT^max = 10. These parameters give converged re- 
sults. If the error from the PP approximation were neg- 
ligible, all EOS curves would superimpose on each other. 
Instead, Fig.[3]shows that the EOS differ significantly. As 
reported in Tab. ITU the transition pressure seems to con- 
verge with respect to the Troullier-Martins core radius as 
it gets small. Its main effect is to shifts the relative posi- 
tion between the diamond and /3-tin EOS branches. The 
EOS from PP's calculations are however different from 
the all-electron one, even for the smallest r c . Our re- 
sults clearly show how the prediction of properties under 
pressure is affected by the PP approximations. The best 
choice in the DFT framework is to use the PAW pseu- 
dopotentials, which give both the transition pressure and 
EOS very close to the all-electron results. 

At the QMC level additional errors may came from 
the lack of a consistent method to generate PP's from the 
corresponding correlated QMC calculation for an isolated 
atom£2. A direct evaluation of the core- valence interac- 
tion was attempted in Refll8l. 

In Tab. |m] we report our VMC and LRDMC results 
using Hartree-Fock, energy adjusted and LDA generated 
pseudopotentials. The transition pressures reported in 



DFT method 


Pi (GPa) 


LDA LAPW 


7.12 


LDA PW PAW 


7.21 


LDA PW Troullier-Martins r c = 1.67 a 


7.44 


LDA PW Troullier-Martins r c = 2.2 an 


7.65 


LDA PW Troullier-Martins r c = 2.6 an 


8.27 



TABLE II: DFT transition pressures for different pseu- 
dopotentials in PWSCF calculations and all-electron Wien2K 
LAPW calculations done with the LDA functional. Note the 
convergence of the Troullier-Martins pseudopotential with re- 
spect to the core radius r c . By using soft pseudopotentials, 
the error could be of 1 GPa on the final transition pressure. 
The best DFT pseudopotential is the PAW one, with a fi- 
nal transition pressure within 0.1 GPa from the all-electron 
result. 
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FIG. 3: EOS obtained by DFT LDA calculations with 
PWSCF and LAPW methods. The PP's in the PW formal- 
ism have been generated by the Troullier-Martins scheme with 
the cutoff radius r c reported in the figure key. The zero of 
the energy has been chosen to be the minimum of the /3-tin 
EOS, corresponding to the volume V — 15 A 3 /Si. This choice 
helps the comparison among the EOS curves. Note that the 
PWSCF calculations done with the PAW pseudopotential are 
on top of the all-electron LAPW points. 
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pseudopotential 


VMC 


LRDMC 


HF Trail-Needs 
energy adjusted Burkatzki-Filippi-Dolg 
LDA Troullier-Martins r c = 1.67 a 
LDA Troullier-Martins r c = 2.1 a 


15.48(6) 
15.80(6) 
15.39(18) 
15.08(18) 


16.65(15) 
16.50(12) 
16.83(19) 
16.11(19) 



TABLE III: Transition pressures for different pseudopoten- 
tials in VMC and LRDMC calculations. Calculations are done 
at volume per atom of 13.08A 3 and 19.95A 3 for the diamond 
and /3-tin phase respectively, and the Maxwell construction 
done assuming for the EOS the same curvature as in the 
Trail-Needs case, fully resolved with respect to its volume 
dependence. r c is the cutoff radius in the pseudopotential 
generation scheme. The factor used to convert energy differ- 
ences into pressures is 26.812 GPa/eV. 

Tab. lIIII are evaluated by computing the energy of the dia- 
mond and /3-tin phase at the volume per atom of 13.08A 3 
and 20.69 A 3 respectively, and assuming that the curva- 
ture of the corresponding EOS is the same as the one 
computed for the Trail-Needs pseudopotentials. 

Results in Fig. |3] and Tab. Mil clearly shows the inter- 
action between core and valence electrons cannot be ap- 
proximated by a rigid shift in energy ( as usually assumed 
estimating the effect of PP on the transition pressure) . 

Core-valence interactions accounts for a correction of 
1.2 ± 0.6 GPa and were calculated in Ref. In order 
to improve the accuracy of this correction an all electron 
calculation is required. At present this is almost impossi- 
ble within QMC, and therefore an uncertainty of at least 
1 GPa coming directly from the PP approximation is un- 
avoidable in our QMC findings. 



B. Phonons and temperature effects 

The inclusion of finite temperature and zero point mo- 
tion effects is crucial for a direct comparison of our results 
with finite temperature experiments (usually performed 
at room temperature). Both experiments^ and theory^ 
indicate that temperature corrections induce a positive 
shift to the critical equilibrium line. On the other hand, 
a further shift in pressure is induced by the inclusion of 
zero point motion effects. Phonon dispersions are in fact 
different in the two phases and zero point motion effects 
do not compensate. 

Finite temperature effects and zero point motion ener- 
gies are included in our estimate of the transition pres- 
sure, by calculating phonon dispersions for the /3-tin and 
diamond phase with the QE package. In order to study 
the convergence of our phonon calculations with the k- 
point mesh we performed PWSCF runs with large plane 
wave cutoff (up to 100 Ry), accurate k-point sampling 
(up to 1620 inequivalent k-points for the /?— tin phase 
and 480 inequivalent k-points for the diamond phase) , as 
well as a small value of the Gaussian broadening (width 
of 0.0001 Ry). We have used the PAW data set as in the 
previous section for the LDA calculation, and an ultra- 



soft (US) PP generated with similar parameters and core 
radii for the PBE calculation. 

Following Ref. 3, we compute the harmonic correction 
AF to the free energy per atom, by using the phonon 
density of states D(lu) available in QE after Fourier in- 
terpolation of the phonon bands, (namely by using mat- 
dyn.x): 

oo 

AF = kT J duD{uj) In (2 sinh (^f ) ) (4) 



where h = 1 is assumed. In using the above expression, 
one has to take into account that the total phonon density 
of states per atom is obviously normalized to 

00 

J dujD(uj) = 3 (5) 


as there are three phonon modes per atom in the thermo- 
dynamic limit. Integrations were changed to summations 
over a uniform mesh with high resolution (lcrra -1 ), and 
the original density of states was appropriately scaled to 
fulfill Eq. El 

Free energy corrections (see Tab. IIV[) are then added 
to a total energy zero temperature calculations. In this 
way our calculation of the transition pressure, estimated 
by the Maxwell construction of the free energy curves, 
is essentially free of systematic errors within the chosen 
DFT functional, as long as phonon anharmonic effects 
can be neglected in the low temperature regime. This is 
a reasonable assumption below the melting temperature 
occurring at about ~ lOOOif. 

To test the impact of the exchange and correlation 
functional on quantum correction estimates, we have per- 
formed calculations with both the PBE and LDA func- 
tionals. As shown in Fig.Hl although different functionals 
provide different transition pressures, the corrections to 
the bare values are very similar and consistent within 
0.2 GPa and in fair agreement with experimental results. 
The results demonstrate that phonons are rather well 
described within DFT and these corrections are very re- 
liable at least before the melting point. Our results do 
not agree with a previous worki^ on this subject, where 
the zero temperature quantum corrections were under- 
estimated by about a factor two, and finite temperature 
corrections were larger by about a factor three. PP used 
111 RefH is no more available and we were not able to 
reproduce the quoted results. Presently the reason for 
this discrepancy is not clear. 

With new available PP's, our temperature corrections 
are very well converged, and appear in reasonable agree- 
ment with recent experimental data (see Fig.[3J. 

By using DFT-PBE free energy corrections to VMC 
total energies, we obtain the corrected VMC curve re- 
ported in Fig. [51 and the corrected transition pressures 
in Tab. [Vj As one can note, the zero point energy for 
the diamond phase is larger than the one for the /3-tin 
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V (A 3 /Si) OK 100K 300K 500K 700K lOOOK 

Diamond phase 

15.882 0.077/0.073 0.073/0.069 0.042/0.039 -0.012/-0.016 -0.084/-0.090 -0.218/-0.226 

17.169 0.069/0.068 0.067/0.066 0.040/0.041 -0.012/-0.011 -0.083/-0.081 -0.215/-0.212 

18.524 0.066/0.064 0.064/0.063 0.039/0.038 -0.012/-0.014 -0.082/-0.084 -0.213/-0.216 

19.228 0.064/0.062 0.062/0.061 0.038/0.036 -0.014/-0.016 -0.084/-0.087 -0.216 -0.220 

19.949 0.062/0.060 0.060/0.059 0.036/0.034 -0.017/-0.019 -0.088/-0.091 -0.221/-0.225 

20.687 0.060/0.058 0.058/0.057 0.033/0.031 -0.020/-0.022 -0.092/-0.095 -0.226/-0.231 

21.444 0.058/0.056 0.056/0.055 0.031/0.029 -0.023/-0.026 -0.096/-0.101 -0.233/-0.238 

23.013 0.054/0.052 0.052/0.051 0.025/0.023 -0.032/0.034 -0.108/0.112 -0.249/0.254 

24.656 0.050/0.048 0.048/0.047 0.019/0.017 -0.040/0.044 -0.120/0.124 -0.265/0.271 

/?— Sn phase 

13.081 0.055/0.053 0.053/0.051 0.025/0.022 -0.033/0.036 -0.110/0.114 -0.251/0.257 

14.004 0.049/0.047 0.048/0.046 0.017/0.014 -0.044/0.048 -0.125/0.130 -0.272/0.280 

15.000 0.044/0.042 0.042/0.040 0.008/0.005 -0.058/0.063 -0.144/0.150 -0.299/0.308 

15.978 0.039/0.037 0.036/0.035 -0.002/0.005 -0.073/0.078 -0.164/0.172 -0.328/0.339 

17.031 0.034/0.032 0.030/0.028 -0.014/0.019 -0.093/0.100 -0.192/0.202 -0.368/0.382 

18.984 0.026/0.025 0.020/0.019 -0.037/0.040 -0.129/0.135 -0.242/0.250 -0.438/0.450 

TABLE IV: Finite temperature corrections to the free energy (eV/atom) within PBE/LDA DFT theory, as explained in the 
text. 
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FIG. 4: Effect of harmonic quantum corrections on the tran- 
sition pressure in Silicon as a function of temperature. 





VMC p t 


LRDMC p t 


Born-Oppenheimer 
zero point motion 
T = 300 K 


15.48(6) 
14.84(6) 
14.53(6) 


16.65(15) 
15.71(14) 
15.74(17) 



TABLE V: Transition pressures p t in VMC and LRDMC cal- 
culations with quantum and temperature effects. The values 
are obtained by a Maxwell construction on a cubic fitting of 
the EOS, corrected at each volume by the DFT-PBE esti- 
mates. 



FIG. 5: The corrected equation of state once the zero point 
energy and the temperature effects are added (estimated at 
the DFT-PBE level). The reference is the VMC calculation 
reported in Fig. [7] and here the zero of energy is taken as the 
VMC value at V = 13.08 A 3 /Si. 

one estimated in Ref.[3 At the LRDMC level, we obtain 
roughly the same total correction (0.91 GPa), although 
it is more difficult to discriminate the temperature effect, 
as the statistical error is larger (see Tab. |V|) . 

C. Kleinman-Bylander approximation 



phase by ~ 0.2 eV, which decreases the transition pres- 
sure by 0.65 GPa. The finite temperature correction is 
negative, and its absolute value is larger for the /3-tin 
phase. Therefore, a further reduction of 0.30 GPa is ob- 
tained at 300K, which implies that the total correction at 
room temperature is 0.95 GPa, a smaller value than the 



The matrix elements of the non-local PP can be eval- 
uated either by direct numerical Gauss-Hermite (GH) 
integration over the polar coordinates or by using the 
Kleinman-Bylander (KB) approximatio n 4 ^ 46 . The KB 
approach is a rather general concept, and it is applied in 
DFT calculations to make the calculation of the PP op- 
erator more efficient. In particular, in the plane- wave 
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formalism the generated PP is conveniently expressed 
in the plane-wave basis set. On the other hand, in 
the QMC framework, one usually works in the coordi- 
nate representation where electron positions and spins 
are given, and this makes the KB construction hard to 
implement numerically. Consequently, the pseudopoten- 
tials used in QMC are usually generated in the so called 
semilocal form (local + non local part), and computed 
by performing a random integration over their angular 
component a 30 ' 47 i 48 . However, in previous QMC calcula- 
tions of the Si diamond-to-/3-tin transitio n 13 ' 18 the deter- 
minantal part has been generated from plane-wave DFT 
calculations, where the KB approximation was used 4 ^ to 
represent a pseudopotential originally written in a semilo- 
cal form. This procedure could lead in principle to a 
poorer form of the variational wave function in the prox- 
imity of the core, where the non-local PP is mostly local- 
ized, with an impact on its nodal structure, and therefore 
a larger FN error. 

To investigate the impact of KB approach, we com- 
pared the total energy, the kinetic energy, and the non- 
local term of the PP, obtained from a single M-point 
DFT calculation with and without the KB approxima- 
tion. The results are shown in Tab. IVII for a system 
of 8 Si in the /3-tin phase with c/a = 0.54 and vol- 
ume 13.08lA 3 / Si, and in the diamond phase with volume 
20.036A 3 /Si. Calculations within the KB approximation 
are done using the PWscf DFT implementation^. We 
use the Qbox codeS 8 - for performing plane-wave calcula- 
tions with Gauss-Hermite integration. An energy cut-off 
of 100 Ry was used in all the plane- wave calculations. 
The exchange and correlation energy was described in 
the LDA by the Perdew-Zunger functional^. The ener- 
gies reported in Tab. IVII and Tab. IVIII clearly show that 
the use of KB approximation causes an error of 0.3 eV/Si 
when evaluating the contribution from the non-local term 
of the PP. This error cancels out in the energy difference 
between the two phases, leaving the DFT results unbi- 
ased by this kind of approximation. In principle, our 
analysis cannot exclude that the nodal structure of the 
DFT generated wave function is unaffected close to the 
core, since there is a significant KB error in the PP con- 
tribution of the total energy. In practice, the close agree- 
ment between our LRDMC results, unaffected by the KB 
approximation, and the ones in Ref. Il3l where this ap- 
proximation has been used, shows that the cancellation 
of KB errors applies also in QMC calculations and leads 
to unbiased energy differences. 

Our results are unaffected by the KB approximation 
by construction, as the DFT code implemented in the 
TurboRVB package uses the exact GH integration of the 
semilocal pseudopotentials. Moreover, the DFT results 
shown in Tab. IVII and Tab. IVIII are a further validation 
of the convergence of our periodic Gaussian basis set. 
Indeed, we found a perfect agreement between the en- 
ergies obtained with plane-wave calculations (performed 
by means of the Qbox code2£ and without the KB ap- 
proximation) and the ones obtained with our Gaussian 
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FIG. 6: Energy vs c/a in the /?— tin phase at fixed volume 
l°ok 3 /Si for various sizes. FS corrections from the KZK ex- 
trapolation are included and energies are averaged over the 
two most symmetric points (r and M) of the Brillouin zone. 
In the largest size the Jastrow factor was parameterized at 
large distance according to Eq.@ and the error of the re- 
stricted variational freedom (~ —0.0082 eV/atom) was esti- 
mated with the same calculation for 64 Si atoms. 



basis set implementation in the TurboRVB packagi 



,29 



IV. RESULTS 

In this section we present our QMC results for the 
transition pressure, performed with the Trail-Needs pseu- 
dopotentials. We first validate the quality of our vari- 
ational wave function in both the /3-tin and diamond 
phases. In the /3-tin phase, we compute the c/a ratio in 
the proximity of the critical pressure (at volume V = 15 
A 3 /atom), where this lattice parameter is experimentally 
known. We evaluate also the bulk modulus and other 
structural properties from the fit of the equation of state 
E = E(V), and compare them with former experimental 
and theoretical values. Finally, we compute the transi- 
tion pressure by performing the Maxwell construction on 
the EOS of both phases. 

In Fig. [6] we report the VMC energies as a function of 
the c/a ratio for the /3-tin phase at different supercell size 
(64, 96, 256 Si atoms) and fixed volume V — 15 A 3 /atom. 
All the energies include FS corrections within the KZK 
scheme, as reported in Subsec fllBI The results in Fig. [6] 
are a first test of the quality of our WF for describing 
the /3-tin phase. With the Jastrow optimization we re- 
produce rather well the experimental value c/a — 0.554. 
Therefore for all the following VMC and LRDMC calcu- 
lations we fix the c/a ratio to the value 0.55. 

In Fig. [7] we reported the VMC results for the energy 
E(V) as a function of the volume V. All the energies in 
the Figure include FS corrections using the KZK scheme. 
Although FS effects are more pronounced in the metallic 
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Quantity 


KB 


GH 


present n — 24 


present n = 10 


Total energy 


-102.85692 


-102.84036 


-102.83975 


-102.83737 


Kinetic energy 


50.95461 


50.83712 


50.83563 


50.83195 


Non-local pseudo 


22.13269 


22.45092 


22.45256 


22.45398 



TABLE VI: Comparison of LDA-DFT energies (eV/Si) obtained within the KB approximation and with the GH numerical 
integration over the polar coordinates of non-local PP matrix elements. We report the total energy, the kinetic energy and the 
non-local contribution of the PP for a system of 8 Si with anti periodic boundary conditions (M— point) in the /3-tin phase 
with c/a — 0.54 and volume 13.08lA 3 /5i and Trail-Needs pseudopotentials. Calculations within the KB approximation are 
performed with PWscf code. Plane- wave calculations with a GH numerical integration are done using Qbox code. We report 
also the energies obtained for the same system using the Gaussian basis set implementation of the LDA-DFT method coded in 
the TurboRVB package (referred to as "present" in the table), n indicates the number of Gaussians for single-particle orbitals 
with the s symmetry. 



Quantity 


KB 


GH 


present n — 24 


present n = 10 


Total energy 


-105.32986 


-105.31911 


-105.31880 


-105.31302 


Kinetic energy 


44.0823 


43.96216 


43.96152 


43.95454 


Non local pseudo 


21.43672 


21.79496 


21.79506 


21.79468 



TABLE VII: Same as in Tab. I VII but for a system of 8 Si with anti periodic boundary conditions in the diamond phase at 
volume 20.036A 3 /Si. 



phase, the results for the /3-tin phase clearly show that 
the FS errors are under control. In fact we find that the 
energies for the 256 atoms supercell fall on the top of 
the data for the 64 atoms calculations. LRDMC ener- 
gies are shown in Fig. [8] for 64 atoms supercell. VMC 
and LRDMC EOS are fitted using a cubic polynomial 
function. The results for bulk properties of the diamond 
phase, reported in Tab. I Villi are in very good agreement 
with the experimental values, whereas the ones for the 
/3— tin structure compares well with previous QMC data. 
Note also that there exists a sizable zero point motion 
correction to both the equilibrium volume and the bulk 
modulus, that has not been taken into account so far in 
previous works. We have estimated these corrections by 
adding to the EOS the zero temperature quantum correc- 
tions evaluated within the harmonic approximation and 
the PBE functional, as explained in Subsec. IIII B[ 

The critical pressure of the diamond to /3-tin transition 
is reported in Tab. IIXI VMC calculations give a raw pt 
of 15.48(6) GPa, LRDMC data give 16.65(15) GPa. 

The inclusion of zero point motion, finite temper- 
ature and core-valence contributions bring the transi- 
tion pressure to 13.33(70) GPa (VMC) and 14.50(70) 
GPa (LRDMC), for a total final shift of -2.15 GPa 
(Tab. lIXp . Zero point motion and thermal corrections 
at 300 K amount to —0.65 GPa and to —0.3 GPa respec- 
tively (they are estimated performing a PBE phonons 
calculations, as explained in Subsec. IIII Bl) . The core- 
valence interaction contribution is — 1.20 ± 0.6 GPa from 
Ref. [l8|, the same value has been used in Refs. Il3lll8ll20l . 
Contributions beyond frozen core approximation arc 
already included at the DFT level trough non-linear 
core corrections in the PP. We observe that previous 
calculation a 18 , 2013 consider a total correction to the raw 
data of —2.5 GPa, because of the different value of the 
zero point energy and finite temperature effects. All pre- 
vious calculations should be increased by 0.35 GPa, for 
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FIG. 7: Equation of state and the metal-insulator transition 
pressure in bulk Silicon obtained by the VMC technique after 
optimizing a Jastrow factor in a localized basis set containing 
2s and 2p Gaussian orbitals per Si. The energies are corrected 
using the KZK correction scheme.— In this way finite size 
effects of the transition pressure appears to be small. The 
lattice value c/a = 0.55 is used for the /3— tin phase. The 
data plotted here are not yet corrected for the quantum and 
temperature effects. See Tab. IIXI for the final pressure which 
includes zero point motion and thermal contributions. 



accounting this difference. . 

Both our corrected VMC and LRDMC values are 
above the experimental range, and remarkably the more 
accurate LRDMC method leads to a transition pressure 
Pt larger than the VMC result. 

To further support the accuracy of our LRDMC calcu- 
lations, we have systematically optimized the molecular 
orbitals by minimizing the VMC energy in presence of 
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The same as in Fig. [7] but for the LRDMC calcula- 



a Jastrow factor, starting from the LDA orbitals, for a 
system of 8 Silicon atoms (see TablXl). The optimization 
of the molecular orbitals allows us to assert the impact 
of the wave function on the final LRDMC results. 

We have found that the LDA orbitals are a quite accu- 
rate starting point, but for the diamond at equilibrium 
geometry the total VMC energy decreases slightly more 
(~ O.OleV/atom) than the one for the (3— tin at high 
pressure, implying an increase in the transition pressure 
of about 0.4 GPa at the VMC level. That is again consis- 
tent with the more accurate LRDMC calculation, which 
is much less affected by the optimization of the wave func- 
tion as clearly shown in Tab. |Xj Altogether these results 
point even more clearly in the direction of a larger zero 
temperature transition pressure. 



V. CONCLUSIONS 

We have performed DFT and QMC calculations for 
the diamond-to-/3— tin transition in Silicon. At the QMC 
level we have proven that it is possible to accurately and 
efficiently describe correlation effects across the metal- 
insulator transition by applying a relatively simple Jas- 
trow factor to a DFT generated Slater determinant. 

We have shown that the estimation of the PP effect 
is one of the most delicate issues in any QMC calcula- 
tion, and represents so far the most significant source of 
systematic error in the transition pressure, amounting to 
about lGPa, as described in Subsec. IIII Al The core- 
valence correlation correction strongly depends on the 
PP, and is almost unpredictable without doing the cor- 
responding all-electron calculation. Moreover, the only 
known QMC estimate for this correction (— 1.2GPa), was 
done in Ref. [H for the Trail- Needs PP, and is affected by 
a large statistical error (~ 0.6 GPa). 

In the present work we have also found a signifi- 



diamond phase 


V cq (A 3 /Si) 


Ecohcsivc (eV/Si) 


B (GPa) 


LDA 


19.77 


5.29 


95.75 


PBE 


20.42 


4.62 


89.4 


VMC 


20.124±0.036° 


4.6003±0.0015 i ' 


102.3±1.4 C 


LRDMC 


20.33±0.1" 


4.6650±0.003'' 


95.78±3 C 


DMC (Ref. .18) 


20.21± 0.03 a 


4.62 ± 0.01 b 


101.4±10 c 


DMC (Ref. .13) 


20.08±0.05 a 


- 


98±7 c 


Exp. 


20.0 <i 


4.62±0.08 e 


99 d 


/3— Sn phase 


V cq (A 3 /Si) 


E co hcsivc (eV/Si) 


B (GPa) 


LDA 


14.92 


5.10 


115.4 


PBE 


15.45 


4.29 


110.7 


VMC 


15.25± 0.05 * 


4.186 ± 0.0013 9 


119.7 ± 3.5 h 


LRDMC 


15.34 ± 0.16 / 


4.211 ± 0.0024 9 


111.3 ± 8.3 h 


DMC (Ref. .13) 


15.31±0.2 / 




98.6±12 h 



TABLE VIII: Comparison of present numerical results with 
previous QMC data and available experiments for the equi- 
librium properties of Silicon in the diamond structure: equi- 
librium volume (V eq ), cohesive energy (E co hcsivc), and bulk 
modulus (B) are reported. They were estimated by a cubic 
interpolation of the E(V) points. The DFT pseudopoten- 
tials used here includes non linear core polarization terms 
and scalar relativistic corrections and are in consistent agree- 
ment with corresponding all-electron calculations. For the 
DFT methods the cohesive energy is given by the magnetic 
solution of the single atom, whereas quantum corrections are 
included according to Tab. (|IV |). All corrections to QMC data 
are estimated within the DFT-PBE functional. We have also 
checked that standard LDA functional provides similar cor- 
rections. 

"Corrected by 0.lA 3 /atom to take into account the zero point 
motion. 

'Corrected by -0.06eV/atom to take into account the zero point 
motion. 

c Corrected by -1.6GPa to take into account the zero point motion. 
VMC and LRDMC are further corrected by -10 GPa to take into 
account non cubic terms in the interpolation within the equilibrium 
volume range 13 o 19 (16 <r> 25) A 3 /atom in the /3-tin (diamond) 
phase. 

d Taken from Ref. pi] 

e Taken from Ref. [53 

^Corrected by 0.11 A 3 /atom to take into account the zero point 
motion. 

9 Corrected by -0.04eV/atom to take into account the zero point 
motion. 

^Corrected by -2.8GPa to take into account the zero point motion. 
VMC and LRDMC are further corrected by -5.6Gpa analogously 
to the diamond phase. 



cant reduction of the phonon correction to the transi- 
tion pressure- from the quoted —1.3 GPa 14 to our fully 
converged LDA and PBE value of -0.95 GPa at 300K. 
This shifts up all previous QMC pt estimates by ~ 0.35 
GPa. Thus, we arrive at the conclusion that all Monte 
Carlo findings published so far predict a transition pres- 
sure that is significantly larger than the one observed at 
room temperature, from 1/2 GPa3£ to 4 GPa 1 ^ above the 
experimental upper edge of 12.5 GPa. 

Our VMC technique is in agreement with the value of 
the transition pressure recently reported by the AFQMC 
technique 20 . A clear increase of the transition pressure 
by about lGPa is obtained when the accuracy of the cal- 
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method 


raw (GPa) 


corrected (GPa) (T = 300 K) 


LDA 


7.21 


6.34 


PBE 


9.87 


8.99 


VMC 


15.48±0.06 


13.33±1.0 


LRDMC 


16.65±0.15 


14.50±1.0 


DMC (Ref. 18) 


19.0 ± 0.5 


16.5±0.5 


DMC (Ref. 13j 


16.5±1.0 


14.0 ± 1.0 


AFQMC (Ref. 20) 


15.1 ± 0.3 


12.6 ± 0.3 


Exp. 


10.0 - 12.5 


10.0 - 12.5 



whereas the DFT Slater determinant was obtained with 
the PBE functional and the KB approximation for the 
pseudopotential. This agreement suggests also that 
the DMC/LRDMC technique is weakly dependent on the 
functional used to generate the DFT orbitals and should 
be considered rather accurate for a given pseudopoten- 
tial. 

To summarize, our best estimate of the transition pres- 
sure at room temperature is 14.5 GPa, with an uncer- 
tainty of 1 GPa coming mainly from the pseudopoten- 
tial approximation, as our finite-size extrapolation error 
is definitely smaller. The discrepancy with respect to 
the experimental values leads us to conclude that further 
work is necessary to determine the phase boundary of 
the metal-insulator transition in Silicon. On one hand, 
from the experimental point of view one should verify 
whether, by removing the stress anisotropy in the exper- 
imental environment, the transition pressure can signifi- 
cantly increase and get closer to the QMC prediction, as 
suggested in Ref.Qjl On the other hand, in QMC calcula- 
tions it should be worth defining consistent pseudopoten- 
tials, since we have seen that they can significantly affect 
the EOS at large pressure. So far there is in fact no accu- 
rate method to estimate the systematic error related to 
the pseudopotential approximation, since an all-electron 
calculation of bulk silicon is basically prohibitive within 
the QMC method. A first attempt along these lines has 
been done in Ref. [H. At this stage of development the 
construction of pseudopotentials is quite unsatisfactory 
for high accuracy QMC calculations, since the pseudopo- 
tentials are usually determined with different and less 
accurate techniques, as Hartree-Fock or LDA. Despite 
the recent progress in the use of pseudopotentials within 
DMC 3 - ' 40 ' 54 , the implementation of the pseudopotential 
approximation in the many-body QMC framework is not 
as mature as in the DFT, where a remarkable progress 
was made only after several years of experience with the 
so called PAW method 4 ^. Thus, in QMC we believe 
there is room for a significant improvement to be real- 
ized in the next years. 



Appendix A: Gaussian periodic basis set 

We use a localized Gaussian basis set on a box of of 
lengths L x , L y , L z , defining the periodic electron- ion dis- 
tance as 



TABLE IX: Zero temperature transition pressure in GPa 
obtained by a cubic interpolation of the EOS. Comparison 
of the present numerical results with available experiments 
and previous theoretical data. The corrected numerical QMC 
data are obtained after including the zero point motion, fi- 
nite temperature, and core-valence contributions, which are 
not present in the raw data, as explained in the text. DFT 
corrections include only the zero point motion and finite tem- 
perature effects as they are performed with the non-linear core 
correction in the pseudopotentials. In particular, the LDA re- 
sults have been obtained by using a PAW PP, while the PBE 
values come from an US PP. For the core-valence contribu- 
tions we have used the published DMC estimate in Refill, 
that unfortunately is affected by a very large error bar, as 
discussed in the text. The raw data and the corrections ap- 
plied by other authors have been taken from the referenced 
papers. 

culation is improved by the LRDMC scheme. Although 
LRDMC total energies are more accurate, it is in prin- 
ciple possible that this technique could lead to results 
worse than the VMC ones in the estimation of the equa- 
tion of state. However we believe that this is quite un- 
likely, because DMC as well as LRDMC always improve 
the accuracy of the physical estimates, as long as they are 
applied to a reasonably good variational wave function. 

As a further independent check that a more accurate 
transition pressure is larger than the VMC estimate, we 
have optimized the Slater determinant in presence of 
the Jastrow factor for a small number of atoms, and 
found a consistent increase in the transition pressure. Al- 
though we cannot estimate more accurately the nodal 
error -namely the exact ground state result for given 
pseudopotential- it looks plausible that by fully optimiz- 
ing the wave function this error should decrease. There- 
fore, our orbital optimization should give at least the 
trend of the correction to the approximate VMC result. 

On the other hand, our LRDMC result is very close 
to the recent DMC calculation performed by Hennig 
et al^. where the same pseudopotential was used, 



= V(v si ° - x »))' + (v si ° (£<» - >»»)) S + - z "jf < A1 > 

I 

where Xi,yi,Zi indicate the Cartesian components of the corresponding ion ones Rj with j — 1,---Na- The 
the electron coordinates r j for i = 1, • • • N and Xj , Yj , Zj 
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system 


DFT+J/LRDMC 


FOPT/LRDMC 


diamond V= 19.949 A 3 /Si PBC 


-106.0120(12) /-106.3064(27) 


-106.0493(10) /-106.3204(24) 


/3-tin V= 13.081A 3 /Si PBC 


-106.5557(11)/-106.8765(28) 


-106.5890(10)/-106.8871(32) 


/3-tiri V= 13.08lA7Si APBC 


-103.7359(13) /-104.0720(41) 


-103.7617(9)/-104.0818(36) 



TABLE X: Variational Monte Carlo/LRDMC energy per atom (eV) for a system of 8 Si with periodic (PBC) and antiperiodic 
(APBC) boundary conditions obtained with the same basis: 8s6p4d for the Slater determinant, and 2s2p/lslp for the Jastrow 
factor. In the DFT+J case only the Jastrow factor (with no restriction to the exponents of the Gaussians) is optimized, while 
the determinantal part is the output of an LDA calculation in the 8s6p4d basis. Conversely, in the fully optimized (FOPT) 
case the molecular orbitals are optimized together with the Jastrow factor, while keeping fixed the exponents of the 8s6p4d 
Gaussians to the even tempered values discussed in the text (see Subsec. [II A[) . 



angular part of the Gaussian basis can be defined in strict 
analogy with the conventional scheme for open systems. 
They are obtained by multiplying the overall Gaussian 
exp(— Zrfj) by appropriate polynomials of sin(-^-r^ 7 ) 

and cos(-^-r-j) , where r^j = r£* — R^ 1 and /i = x,y,z 
labels the three Cartesian components. Strictly speaking 
in the periodic case there is an arbitrariness in defining 
this polynomials because the multiplication of a polyno- 
mial by any even cos-power defines an allowed element 
of the atomic basis satisfying the periodic or antiperiodic 
boundary conditions in each direction. 

In order to define the basis in an accurate and conve- 
nient way, we have chosen the unique polynomials that 
contain the minimum possible cosine powers. For in- 
stance the angular part Yi = 2,mf 2 of the <i-wave orbital 
is defined by: 

for m = : (A2) 

^( 3 (^ Sin( £ r3) ) 2 - r2 ) 

for m = ±1 : (A3) 
sin^r 2 ) cos(£r 2 )^ sin(£r 3 ) cos(£r 3 ) 
sin^r 1 ) cos^r 1 )-^ s i n (-^-r 3 ) cos(-^r 3 ) 

for m = ±2 : (A4) 
sin^r 1 ) cos(f r r 1 )^ sin(£r 2 ) cos(£r 2 ) 

#((^^(fr rl )) 2 -(v^(^ 
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